#import arcgisscripting, sys, os, glob
import time as timemod
#import math
from numpy import *
#from pylab import *
from ArcGISinterface import *
from pyproj import Proj
#from scipy.io import write_array
from LatLongUTMconversion import LLtoUTM

#class GPSdata(IsDescription):
#   profile = Int64Col()   
#   trace = Int64Col()     
#   northing = FloatCol()    
#   easting = FloatCol()
#   date = StringCol(16)
#   time = StringCol(16)

def wxOpenFile(dir=os.getcwd()):
    import wx
    
    extfilter = "All Files (*.*)|*.*|C-band GPS Files (*.gps)|*.gps"
    app = wx.PySimpleApp()

    dlg = wx.FileDialog(None, message='Select files ...', defaultDir=dir, wildcard=extfilter, style=wx.FD_MULTIPLE)  
    if dlg.ShowModal() == wx.ID_OK:
        pathlist = dlg.GetPaths()
        filelist = []
        for path in pathlist:
            filelist.append(os.path.basename(str(path)))
        filelist.sort()
        os.chdir(os.path.dirname(str(pathlist[0])))
        return filelist
    else:
        return None
    dlg.Destroy()


def cbandgps(filename):
    ''' Reads '.gps' files from C-band GPR '''
    fid=open(filename,'r')
    text=fid.readlines()
    fid.close()
#    map = Basemap(projection='spstere',boundinglat=0,lon_0=0, lat_0=-71,resolution='l')
    
    project = Proj(proj='stere', ellps='WGS84', lat_0=-90.0, lon_0=0.0, lat_ts=-71.0, x_0=0.0, y_0=0.0)

    N=len(text)
    traceNo=[]
#    Lat=[]
#    Lon=[]
    X = []
    Y = []
    recorded=[]
    profileNo=[]
    pctime=[]
    for n in range(1,N):
        coords=asarray(text[n].split(' '))
        if float(coords[1])!=0.:
            traceNo.append(int(coords[0]))
            x, y = project(coords[2], coords[1])
#            Lat.append(float(coords[1]))
#            Lon.append(float(coords[2]))
            X.append(x)
            Y.append(y)
            recorded.append(timestruct(coords[3], coords[4]))
            pctime.append(float(coords[5]))
          
    return traceNo, Y, X, recorded, pctime


def pppgps(filename):
    ''' reads Precision Point processed .csv files '''
    fid=open(filename,'r')
    text=fid.readlines()
    fid.close()
    N=len(text)/2
    Lat=[]
    Lon=[]
    Hgt=[]
    Hour=[]
    Day=[]
    Year=[]
    for n in range(1,N):
        coords=asarray(text[n].split(','))
        if float(coords[1])!=0.:
            Lat.append(float(coords[0]))
            Lon.append(float(coords[1]))
            Hgt.append(float(coords[2]))
            Hour.append(float(coords[3]))
            Day.append(int(coords[4]))
            Year.append(int(coords[5]))
            
    return Lat, Lon, Hgt, Hour, Day, Year
    
    
def travgps(filename):
    fid=open(filename,'r')
    text=fid.readlines()
    fid.close()
    N=len(text)
    Lat=[]
    Lon=[]
    Name=[]
    Date=[]
    for n in range(1,N):
        coords=asarray(text[n].split(';'))
        Lat.append(float(coords[6]))
        Lon.append(float(coords[7]))
#        Alt.append(float(coords[8]))
        Name.append(coords[2])
        Date.append(coords[4])
            
    return Name, Date, Lat, Lon#, Alt
    

def gpsplot(E,N):
    scatter(E,N)
    xlabel('Easting')
    ylabel('Northing')
    show()

def array2mdb(feat,traceNo, Lat, Lon, recorded, profileNo):
    gp = arcgisscripting.create()
    gp.Overwriteoutput = 1
    
    rows = gp.InsertCursor(feat)
    for n in range(len(Lat)):
        pnt = gp.CreateObject("Point")
        #pnt.id = n+1
        pnt.x = Lon[n] 
        pnt.y = Lat[n] 
        row = rows.NewRow()
        row.shape = pnt
        row.profileNo = profileNo
        row.Lat = pnt.y
        row.Lon = pnt.x
        rec = '%s.%s.%s' % (recorded[n][2], recorded[n][1], recorded[n][0])
        row.recorded = rec
        row.traceNo = traceNo[n] 
        rows.InsertRow(row)
        
def csv2mdb(feat, Lat, Lon, Hgt, Hour, Day, Year):
    gp = arcgisscripting.create()
    gp.Overwriteoutput = 1
    
    rows = gp.InsertCursor(feat)
    for n in range(len(Lat)):
        pnt = gp.CreateObject("Point")
        #pnt.id = n+1
        pnt.x = Lon[n] 
        pnt.y = Lat[n] 
        pnt.z = Hgt[n]
        row = rows.NewRow()
        row.shape = pnt
        row.Lat = pnt.y
        row.Long = pnt.x
        row.Elev = pnt.z 
        row.recorded = converttime(Year[n], Day[n], Hour[n])
        rows.InsertRow(row) 
        
def trav2mdb(feat, Name, Date, Lat, Lon):
    gp = arcgisscripting.create()
    gp.Overwriteoutput = 1
    
    rows = gp.InsertCursor(feat)
    for n in range(len(Lat)):
        pnt = gp.CreateObject("Point")
        #pnt.id = n+1
        pnt.x = Lon[n] 
        pnt.y = Lat[n] 
#        pnt.z = Alt[n]
        row = rows.NewRow()
        row.shape = pnt
        row.Lat = pnt.y
        row.Lon = pnt.x
#        row.Alt = pnt.z
        row.Name = Name[n]
        row.Datum = Date[n] 
        rows.InsertRow(row)   
        

def converttime(yy,dd,hh):
    '''
    input from pppgps
    output time as hh:mm:ss date as 
    '''
    hour=int(hh)
    if hour==24:
        hour=0
        dd=dd+1
    mm=hh%1*60.0
    min=int(mm)
    ss=mm%1*60.0
    sec=int(ss)
    timestr='%i %i %i %i %i' % (yy, dd, hour, min, sec)
    tstruct=timemod.asctime(timemod.strptime(timestr, '%Y %j %H %M %S'))
    return tstruct

def timestruct(date,time):
    date=date.split('/')
    time = time.split(':')
    ts = timemod.strptime(date[2]+' '+date[1]+' '+date[0]+' '+time[0]+' '+time[1]+' '+time[2], '%Y %m %d %H %M %S')
#    ts = (int(date[2]), int(date[1]), int(date[0]), int(time[0]), int(time[1]), int(time[2]), int(0), int(1), int(-1))
    return ts

def writegps(gps, dir, filename): 
    fid=open(filename+'.txt','w')
    for i in range(len(gps)):
        fid.write('%s %f %f\n' % (dir+'_'+filename,gps[i,0],gps[i,1]))
    fid.close()

def interpolateDEGgps(traceNo, Lat, Lon, recorded, totaltraces):
    '''
    Adds n points between two GPS points
    '''       
    DEGgps=gprgps()
    N=len(Lat)
    
    # insert coords for traces 0 to first traces with GPS measurement 
    N1, E1=Lat[0], Lon[0]
    N2, E2=Lat[1], Lon[1]
    dE=E2-E1
    dN=N2-N1
    tdiff=traceNo[0]
    for i in range(tdiff):
        DEGgps.Easting.append(E1+dE/tdiff*i)
        DEGgps.Northing.append(N1+dN/tdiff*i)
        DEGgps.traceNo.append(i)
        
    # insert coords for traces from first to last GPS measurement
    for i in range(N-1):
        N2, E2=Lat[i+1], Lon[i+1]
        tdiff=traceNo[i+1]-traceNo[i]
        dE=E2-E1
        dN=N2-N1
        for k in range(tdiff):
            DEGgps.Easting.append(E1+dE/tdiff*k)
            DEGgps.Northing.append(N1+dN/tdiff*k)
            DEGgps.traceNo.append(traceNo[i]+k)
        E1=E2
        N1=N2
        
    # insert coords for last GPS measurement until last trace
    E1, N1=Lat[-2], Lon[-2]
    E2, N2=Lat[-1], Lon[-1]
    dE=E2-E1
    dN=N2-N1
    tdiff=totaltraces-traceNo[-1]
    for i in range(tdiff):
        DEGgps.Easting.append(E1+dE/tdiff*i)
        DEGgps.Northing.append(N1+dN/tdiff*i)
        DEGgps.traceNo.append(i)
        
    return DEGgps

def interpolateUTMgps(traceNo, Lat, Lon, recorded, totaltraces):
    '''
    Adds n points between two GPS points
    '''       
    UTMgps=gprgps()
    N=len(Lat)
    
    # insert coords for traces 0 to first traces with GPS measurement 
    [UTMZ1, N1, E1]=LLtoUTM(23, Lat[0], Lon[0])
    [UTMZ2, N2, E2]=LLtoUTM(23, Lat[1], Lon[1])
    dE=E2-E1
    dN=N2-N1
    tdiff=traceNo[0]
    for i in range(tdiff):
        UTMgps.UTMzone.append(UTMZ1)
        UTMgps.Easting.append(E1+dE/tdiff*i)
        UTMgps.Northing.append(N1+dN/tdiff*i)
        UTMgps.traceNo.append(i)
        
    # insert coords for traces from first to last GPS measurement
    for i in range(N-1):
        [UTMZ2, N2, E2]=LLtoUTM(23, Lat[i+1], Lon[i+1])
        tdiff=traceNo[i+1]-traceNo[i]
        dE=E2-E1
        dN=N2-N1
        for k in range(tdiff):
            UTMgps.UTMzone.append(UTMZ1)
            UTMgps.Easting.append(E1+dE/tdiff*k)
            UTMgps.Northing.append(N1+dN/tdiff*k)
            UTMgps.traceNo.append(traceNo[i]+k)
        UTMZ1=UTMZ2
        E1=E2
        N1=N2
        
    # insert coords for last GPS measurement until last trace
    [UTMZ1, N1, E1]=LLtoUTM(23, Lat[-2], Lon[-2])
    [UTMZ2, N2, E2]=LLtoUTM(23, Lat[-1], Lon[-1])
    dE=E2-E1
    dN=N2-N1
    tdiff=totaltraces-traceNo[-1]
    for i in range(tdiff):
        UTMgps.UTMzone.append(UTMZ1)
        UTMgps.Easting.append(E1+dE/tdiff*i)
        UTMgps.Northing.append(N1+dN/tdiff*i)
        UTMgps.traceNo.append(i)
        
    return UTMgps


def interpolateDistances(dist, maxTrace):
    '''
    Interpolate distances in between to GPS points.
    Resulting in a distance between each trace in a profile.
    '''
    
    
def filterDuplicates(Lat,Long,Elev,mindist=1.0):
    '''Removes all GPS points that are closer than 'mindist' [m].
    Reduces data and removes stop points in the kinematic GPS data.'''
    N = len(Lat)
    nLat = []
    nLong = []
    nElev = []
    nLat.append(Lat[0])
    nLong.append(Long[0])
    nElev.append(Elev[0])
    duplicates = 0
    for n in range(N-1):
        dist = distance(nLat[-1], nLong[-1], Lat[n+1], Long[n+1])
        if dist > mindist:
#        if Lat[n+1] != nLat[-1] and Long[n+1] != nLong[-1]:
            nLat.append(Lat[n+1])
            nLong.append(Long[n+1])
            nElev.append(Elev[n+1])
#            print nLat[-1], nLong[-1], nElev[-1], dist
        else:
            duplicates = duplicates+1
#    print duplicates        
    return nLat, nLong, nElev

def correctTraceNo(Y, X, recorded, traceEnd, traceInt):
    ''' 
    Determines the trace number that belongs to the measured GPS point
    by dividing the time difference between GPS measurements by the sampling 
    time increment.
    '''
    traceStart = traceInt
    if traceInt == 20:
        samplingtime = 0.25
    else:
        samplingtime = 0.2
    pause_threshold = int(180/samplingtime) # If the time between to measurements is Xer the 3 min a stop is assumed
    traceNo = []
    traceNo.append(traceStart)
    dist = []
    dist.append(0.0)
    distEach = []
    YEach = []
    XEach = []
    # Insert trace distances
    for n in range(len(recorded)-1):
        t1 = timemod.mktime(recorded[n])
        t2 = timemod.mktime(recorded[n+1])
        traceAdd = int((t2-t1)/samplingtime)
        if traceAdd < pause_threshold:
            traceNo.append(traceNo[-1]+traceAdd)
        else:
            traceNo.append(traceNo[-1]+traceStart)
#        dist.append(distance(Y[n], X[n], Y[n+1], X[n+1]))
        dist.append(np.sqrt((X[n+1]-X[n])**2 + (Y[n+1]-Y[n])**2))
        dE = X[n+1]-X[n]
        dN = Y[n+1]-Y[n]
        tdiff = traceNo[-1]-traceNo[-2]
        for m in range(tdiff):
            distEach.append(dist[-1]/tdiff)
            YEach.append(Y[n]+dN/tdiff*m)
            XEach.append(X[n]+dE/tdiff*m)
            
    # Insert the first traces
    dE = X[1]-X[0]
    dN = Y[1]-Y[0]
    for m in range(traceStart):
        distEach.insert(m, (dist[1]/traceStart))
        YEach.insert(m, (Y[0]-dN/traceStart*m))
        XEach.insert(m, (X[0]-dE/traceStart*m))
        
    # Insert the last traces
    dE = X[-2]-X[-1]
    dN = Y[-2]-Y[-1]
    tdiff = int(traceEnd-len(distEach))
    for m in range(tdiff):
        distEach.append(dist[-1]/tdiff)
        YEach.append(Y[-1]-dN/tdiff*m)
        XEach.append(X[-1]-dE/tdiff*m)
        
    distWhole = zeros(len(distEach))
    for i in range(1, len(distEach)):
        distWhole[i] = distWhole[i-1] + distEach[i]
    traceNew = arange(1, traceEnd+1)
    return traceNew, distEach, YEach, XEach, distWhole


def rad2deg(rad):
    """Convert radians to degrees.
    Return float radians."""
    return rad * 180 / math.pi

def deg2rad(deg):
    """Convert degrees to radians.
    Return float degrees"""
    return deg * math.pi / 180

def distance(lat1, lon1, lat2, lon2):
    """Calculate distance between two lat/long pairs.
    Return float distance in meters."""
    lat1 = deg2rad(lat1)
    lon1 = deg2rad(lon1)
    lat2 = deg2rad(lat2)
    lon2 = deg2rad(lon2)
    theta = lon1 - lon2
    dist = math.sin(lat1) * math.sin(lat2) \
         + math.cos(lat1) * math.cos(lat2) * math.cos(theta)
    dist = math.acos(dist)
    dist = rad2deg(dist)
    meters = abs(dist * 60.0 * 1852) # the factor 1852 converts from nautical miles to meters
    return meters


class gprgps:
    def __init__(self):
        self.UTMzone=[]
        self.Easting=[]
        self.Northing=[]
        self.traceNo=[]
        self.recorded=None
        self.Elevation=None
       

#########################################
### Test scripts started in __main__: ###
#########################################

def make_ppp_feature():    
    feat = 'M:\\home\\GIS\\GPS\\PPP_v3.shp'
    dir='M:\\home\\GIS\\GPS\\PPPprocessed\\'

    os.chdir(dir)
    gpsfile=glob.glob('*.csv')
    n=0
    for gpsfile in gpsfile:
        [Lat, Lon, Hgt, Hour, Day, Year]=pppgps(dir+gpsfile)
        [Lat, Lon, Hgt] = filterDuplicates(Lat,Lon,Hgt)
        csv2mdb(feat, Lat, Lon, Hgt, Hour, Day, Year)
        n=n+1; print n
        

def make_gps_feature():
    feat = 'M:\\home\\GIS\\cbandGPS\\cbandGPS_09.shp'
    gpsdir = 'H:\\Data2008_09\\Cband\\processed\\gps_main'
    
    os.chdir(gpsdir)
    gpsfiles = glob.glob('*.gps')
    for file in gpsfiles:
        profileNo = int(file[7:9])
        [traceNo, Lat, Long, recorded, pctime] = cbandgps(file)
        [Lat, Long, Elev] = filterDuplicates(Lat, Long, zeros(len(Lat)), 2.0)
        array2mdb(feat, traceNo, Lat, Long, recorded, profileNo) 
         

        
def test_dist_calc():
    fn = 'H:\\Data2007_08\\CBAND\\071221\\profile10'
    [traceNo, Lat, Lon, recorded, pctime] = cbandgps(fn+'.gps')
    [traceNew, dist] = correctTraceNo(Lat, Lon, recorded)
    
    hist(dist, bins=20)
    figure(2)
    plot(traceNo)
    hold(True)
    plot(traceNew)
    print traceNo[-1], traceNew[-1]
    for n in range(len(dist)-1):
        dt=traceNew[n+1]-traceNew[n]
        print dist[n]/dt
    show()
    
def getProfileDistance():
    import processGPR
    procdir = 'H:\\Data2007_08\\Cband\\processed\\NUS07-5-7'
    os.chdir(procdir)
    filelist = wxOpenFile(procdir)
    filelist.sort()
    print filelist
    traceWhole = asarray([0])
    distWhole = asarray([0.0])
    for file in filelist:
        pro = processGPR.CTD()
        pro.readH5(file[:-3]+'h5')
        [traceNo, Lat, Lon, recorded, pctime] = cbandgps(file)
        [traceNew, distEach, dist] = correctTraceNo(Lat, Lon, recorded, pro.traces)
        traceNew = traceNew + traceWhole[-1]
        dist = dist + distWhole[-1]
        traceWhole = concatenate((traceWhole, traceNew), axis=0)
        distWhole = concatenate((distWhole, dist), axis=0)
    print distEach, len(distEach), traceWhole[-1], distWhole[-1]    
    savez('NUS07-5-7_3gps.npz', traceNo=traceWhole, dist=distWhole)

def getGPSelev():
    
    fn = 'Y:\\GPS\\HMS2007\\071204_071229ed.txt'
    text = loadtxt(fn)
#    Lat = asarray(text[:,4])
#    Long = asarray(text[:,5])
    Elev = asarray(text[:,6])
    x = asarray(text[:,7])
    y = asarray(text[:,8])
    d = zeros((len(x),))
    d[0] = 0.0
    for n in range(len(x)-1):
#        d[n+1] = distance(Lat[n], Long[n], Lat[n+1], Long[n+1])
        d[n+1] = sqrt((x[n+1]-x[n])**2 + (y[n+1]-y[n])**2) + d[n]
    plot(d, Elev)
    print d.max()
    savez('Y:\GPS\HMS2007\\gps.npz', elev=Elev, dist=d)
    show()
    
if __name__ == '__main__':
#    test_dist_calc()
#    make_gps_feature()
#    make_ppp_feature()
#    getProfileDistance()
    getGPSelev()